Investigate if there is a shift towards asexual or sexual reproduction along an altitudinal gradient (altitude range : 1088 - 3212 m.a.s.l.).
Multivariate analyses (ordination) to investigate the data and linear mixed-effects model to describe the change of a ponderate ratio of bulbils with altitude.
Germination of a subset of bulbils and compare the bulbils to flowers ratio of the (clonal) plant derived from the bulbil with that of the parent plant => dead end…
#Height~altitude
#Number of bulbils ~ altitude
#Number of flowers ~ altitude
#"nratioB" ~ altitude
#"reprod" ~ altitude
#"ratiopond" ~ altitude
#"nratioB" ~ reprod
#Ratio of bubils (number)
op<-par(mar=rep(2,4))
coplot(nratioB~altitude|site,data=all)
dotchart(all$nratioB,groups=all$site)
par(op)
#Ponderate ratio
op<-par(mar=rep(2,4))
coplot(ratiopond~altitude|site,data=all)
dotchart(all$ratiopond,groups=all$site)
par(op)
#Ponderate ratio (transformed)
op<-par(mar=rep(2,4))
coplot(ratiopond.tr~altitude|site,data=all)
dotchart(all$ratiopond.tr,groups=all$site)
par(op)
Kendall or Spearman ? In significance testing these two non-parametric coefficients usually produce very similar results and there is no strong reason for preferring one over the other. However, from the calculation point of view, Kendall’s coefficient is usually considered to be the more difficult. On the other hand Kendall’s coefficient does have an advantage over Spearman’s in the respect that its distribution approaches normality more rapidly. Although the two coefficients produce similar results, Spearman’s coefficient tends to be larger than Kendall’s in absolute value (Colwell & Gillett 1982).
#### Fan & Yang (2009) : r=-0.808, P=0.001
#### Bauert (1993) : r=-0.55, P<=0.05
cor.test(all$h,all$altitude, method=c("pearson"))
##
## Pearson's product-moment correlation
##
## data: all$h and all$altitude
## t = -34.3215, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6667833 -0.6106045
## sample estimates:
## cor
## -0.639547
#### Fan & Yang (2009) : r=-0.624, P=0.23
#### Bauert (1993) : rs=-0.75,P<=0.01
cor.test(all$nB,all$altitude, method=c("pearson"))
##
## Pearson's product-moment correlation
##
## data: all$nB and all$altitude
## t = -19.1449, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4592483 -0.3810727
## sample estimates:
## cor
## -0.4209418
#### Fan & Yang (2009) : r=-0.483, P=0.095
#### Bauert (1993) : rs=-0.62, P<=0.05
cor.test(all$nF,all$altitude, method=c("kendall"))
##
## Kendall's rank correlation tau
##
## data: all$nF and all$altitude
## z = -0.7804, p-value = 0.4351
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.01385713
cor.test(all$nF,all$altitude, method=c("spearman"))
##
## Spearman's rank correlation rho
##
## data: all$nF and all$altitude
## S = 839766999, p-value = 0.4488
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.01836023
#### Fan & Yang (2009) : r=-0.617, P=0.025
cor.test((all$nF+all$nB),all$altitude, method=c("pearson"))
##
## Pearson's product-moment correlation
##
## data: (all$nF + all$nB) and all$altitude
## t = -20.9531, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4897864 -0.4142535
## sample estimates:
## cor
## -0.452832
#### Bauert (1993) : rs=0.49, P=0.07
cor.test(all$nratioB,all$altitude, method=c("pearson"))
##
## Pearson's product-moment correlation
##
## data: all$nratioB and all$altitude
## t = -4.8371, df = 1702, p-value = 1.436e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.16303552 -0.06934777
## sample estimates:
## cor
## -0.1164507
##
## Kendall's rank correlation tau
##
## data: all$nF and all$nB
## z = -17.0521, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.3037251
##
## Spearman's rank correlation rho
##
## data: all$nF and all$nB
## S = 1157626778, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4038192
##
## Pearson's product-moment correlation
##
## data: (all$lF + all$lB) and all$h
## t = 51.7256, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7626159 0.7995926
## sample estimates:
## cor
## 0.7817906
##
## Pearson's product-moment correlation
##
## data: sqrt(all$lF + all$lB) and all$h
## t = 51.1701, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7590690 0.7965344
## sample estimates:
## cor
## 0.7784942
##
## Kendall's rank correlation tau
##
## data: all$lF and all$h
## z = 1.7265, p-value = 0.08426
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.03091922
##
## Pearson's product-moment correlation
##
## data: all$lB and all$h
## t = 43.025, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6982419 0.7437884
## sample estimates:
## cor
## 0.7217957
##
## Pearson's product-moment correlation
##
## data: sqrt(all$lB) and all$h
## t = 42.2567, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6915452 0.7379465
## sample estimates:
## cor
## 0.7155342
Theoretical aspects : Legendre & Legendre (2012), Borcard et al. (2011)
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca<-rda(var_env,scale=TRUE))
## Call: rda(X = var_env, scale = TRUE)
##
## Inertia Rank
## Total 10
## Unconstrained 10 10
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 5.9365397 2.0299400 1.0829172 0.6837256 0.1890309 0.0404861 0.0231301
## PC8 PC9 PC10
## 0.0104184 0.0031740 0.0006379
#summary(var_pca)
scores(var_pca,display=c("species"))
## PC1 PC2
## sfroy 1.2495162 0.4450014
## ddeg -1.4170526 -0.4044422
## etpt7 -1.3227407 -0.6978898
## mbal7 1.4061239 -0.3934450
## pday 0.9773010 -1.1247710
## swb 0.8874390 -0.2389113
## avetemp7 -1.4798262 -0.3107048
## cloud7 1.4648185 0.3370283
## prec7 0.8058167 -1.2535949
## srad7 0.1570868 -0.8345368
## attr(,"const")
## [1] 4.864599
par(mfrow=c(1,2))
biplot(var_pca, scaling=1, main="PCA - Distance biplot")
biplot(var_pca,main="PCA - Correlation biplot") #by default scaling 2
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca5<-rda(var_env5,scale=TRUE))
## Call: rda(X = var_env5, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.87589 1.09846 0.77467 0.22187 0.02912
biplot(var_pca5, scaling=1, main="PCA - Distance biplot")
biplot(var_pca5,main="PCA - Correlation biplot") #by default scaling 2
#summary(var_pca5)
scores(var_pca5,display=c("species"))
## PC1 PC2
## sfroy 1.494547 0.4507193
## ddeg -1.540416 -0.6211337
## mbal7 1.727851 0.0936202
## pday 1.350566 -0.7988239
## srad7 0.456630 -1.5621458
## attr(,"const")
## [1] 4.090623
## PC1 PC2 PC3 PC4 PC5
## 2.87588670 1.09845652 0.77466999 0.22186511 0.02912167
## PC1 PC2
## 2.875887 1.098457
## j p
## 1 1 4.00000
## 2 2 9.00000
## 3 3 15.66667
## 4 4 25.66667
## 5 5 45.66667
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca5<-rda(var_env5,scale=TRUE))
## Call: rda(X = var_env5, scale = TRUE)
##
## Inertia Rank
## Total 6
## Unconstrained 6 6
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 3.70583 1.17439 0.81277 0.25510 0.03072 0.02120
biplot(var_pca5, scaling=1, main="PCA - Distance biplot")
biplot(var_pca5,main="PCA - Correlation biplot") #by default scaling 2
#summary(var_pca5)
scores(var_pca5,display=c("species"))
## PC1 PC2
## sfroy 1.4868236 0.3252721
## ddeg -1.5854018 -0.4882835
## mbal7 1.5928166 -0.2031974
## pday 1.1226474 -1.0403298
## srad7 0.3346198 -1.3920851
## altitude 1.6395054 0.4267452
## attr(,"const")
## [1] 4.28139
Escofier & Pagès (2008)
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))
# First import and standardize the data.
# Run PCA and display the biplots :
(all_mes.pca<-rda(all_mes,scale=TRUE))
## Call: rda(X = all_mes, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.75225 1.73346 0.27113 0.17237 0.07079
par(mfrow=c(1,2))
biplot(all_mes.pca, scaling=1, main="PCA - Distance biplot")
biplot(all_mes.pca,main="PCA - Correlation biplot") #by default scaling 2
PCA to detect intercorrelations among traits (St’astná et al. 2012).
# First import and standardize the data.
# Run PCA and display the biplots :
(pop_mes.pca<-rda(pop_mes,scale=TRUE) )
## Call: rda(X = pop_mes, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.92878 1.71005 0.20928 0.10199 0.04991
par(mfrow=c(1,2))
biplot(pop_mes.pca, scaling=1, main="PCA - Distance biplot")
biplot(pop_mes.pca,main="PCA - Correlation biplot") #by default scaling 2
#With Pearson correlation coefficient :
(env.pearson<-cor(scat_env))
## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.8309154 -0.9700733 0.75238565 0.3899248
## sfroy 0.83091541 1.0000000 -0.7617759 0.67438543 0.2874656
## ddeg -0.97007332 -0.7617759 1.0000000 -0.72126970 -0.3604668
## mbal7 0.75238565 0.6743854 -0.7212697 1.00000000 0.8103262
## pday 0.38992478 0.2874656 -0.3604668 0.81032615 1.0000000
## srad7 0.07973121 0.1679037 -0.0395407 0.06109418 0.3251788
## srad7
## altitude 0.07973121
## sfroy 0.16790373
## ddeg -0.03954070
## mbal7 0.06109418
## pday 0.32517881
## srad7 1.00000000
(all.pearson<-cor(scat_all))
## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.8309154 -0.9700733 0.75238565 0.38992478
## sfroy 0.83091541 1.0000000 -0.7617759 0.67438543 0.28746563
## ddeg -0.97007332 -0.7617759 1.0000000 -0.72126970 -0.36046684
## mbal7 0.75238565 0.6743854 -0.7212697 1.00000000 0.81032615
## pday 0.38992478 0.2874656 -0.3604668 0.81032615 1.00000000
## srad7 0.07973121 0.1679037 -0.0395407 0.06109418 0.32517881
## ratiopond 0.45492422 0.3338513 -0.4269430 0.31266283 0.05655636
## srad7 ratiopond
## altitude 0.07973121 0.45492422
## sfroy 0.16790373 0.33385127
## ddeg -0.03954070 -0.42694301
## mbal7 0.06109418 0.31266283
## pday 0.32517881 0.05655636
## srad7 1.00000000 -0.18612483
## ratiopond -0.18612483 1.00000000
#with Kendall coeff. for nratioB
(all.k<-cor(scat_k,method="kendall"))
## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.97960508 -0.84012539 0.56513933 0.2949488
## sfroy 0.97960508 1.00000000 -0.83526865 0.57134088 0.3049593
## ddeg -0.84012539 -0.83526865 1.00000000 -0.55381137 -0.2703157
## mbal7 0.56513933 0.57134088 -0.55381137 1.00000000 0.6728150
## pday 0.29494881 0.30495932 -0.27031572 0.67281501 1.0000000
## srad7 0.02758621 0.03074994 -0.03510972 -0.01887993 0.1484468
## nratioB -0.03076927 -0.04274042 0.01193094 -0.16640444 -0.2207476
## srad7 nratioB
## altitude 0.02758621 -0.03076927
## sfroy 0.03074994 -0.04274042
## ddeg -0.03510972 0.01193094
## mbal7 -0.01887993 -0.16640444
## pday 0.14844676 -0.22074758
## srad7 1.00000000 -0.17394056
## nratioB -0.17394056 1.00000000
#Kendall for
#Plot a matrix of bivariate scatter
source("panelutils.R")
pairs(scat_all,panel=panel.smooth, diag.panel=panel.hist)
pairs(scat_env,panel=panel.smooth, diag.panel=panel.hist)
(Dufour 2009)
pca1<-dudi.pca(betw_env,scann=F,nf=2)
pca1$eig
## [1] 2.87588670 1.09845652 0.77466999 0.22186511 0.02912167
sum(pca1$eig)
## [1] 5
cumsum(pca1$eig)/sum(pca1$eig)
## [1] 0.5751773 0.7948686 0.9498026 0.9941757 1.0000000
#inertie<-cumsum(pca1$eig)/sum(pca1$eig)
#barplot(pca1$eig)
#env. variables per "pop" -> do the bca betw. "site" (not betw. "pop", no sense)
bet1<-bca(pca1,betw$site,scan=FALSE,nf=2) #site as categorical variable
#names(bet1) # to see available bet1$...
bet1$tab
## sfroy ddeg mbal7 pday srad7
## A -0.46031114 0.60126812 0.04698396 0.7841500 0.43422372
## Lie -0.36255505 0.35736913 0.15086431 0.3929474 0.39074027
## Lou -0.12891799 -0.66274369 0.54712982 0.2735975 -1.17880021
## M -0.35557247 0.13641194 -0.22411840 0.2735975 0.90257693
## Sim -0.44686968 0.09785019 -0.65357495 -1.2580599 -0.06856215
## Sio -0.08883799 0.44599064 0.49825820 0.8924490 0.08728575
## V -0.15098294 0.27080858 -0.45721576 -0.6812019 -0.61795459
## ZRo 1.52064623 -1.07153076 0.10019097 -0.5220687 0.24647090
## ZWi 0.28654540 -0.51012378 -0.37512495 -1.0631217 -0.73509981
#bet1$lw : the weights of the groups are the relative frequencies of the categorical variable
#summary(betw$site)/length(betw$site)
bet1$ratio #the between group variance is the inertia of the between PCA
## [1] 0.3395425
#The between inertia is here to 0.3395 i.e. 33.95% of the total inertia is due to the factor.
plot(bet1) #on retrouve le groupe A/Sio/Lie/M
pca2<-dudi.pca(betw_mes,scann=F,nf=2)
pca2$eig
## [1] 2.92877500 1.71004730 0.20927962 0.10199001 0.04990806
sum(pca2$eig)
## [1] 5
cumsum(pca2$eig)/sum(pca2$eig)
## [1] 0.5857550 0.9277645 0.9696204 0.9900184 1.0000000
#inertie<-cumsum(pca2$eig)/sum(pca2$eig)
#barplot(pca2$eig)
#env. variables per "pop" -> do the bca betw. "site" (not betw. "pop", no sense)
bet2<-bca(pca2,betw$site,scan=FALSE,nf=2) #site as categorical variable
#names(bet2) # to see available bet1$...
bet2$tab
## nF nB lF lB h
## A 0.6530908 0.15494952 0.69919429 0.33303322 0.5293246
## Lie 0.1549427 0.52266584 0.21979809 0.23909482 0.5021346
## Lou 0.1056213 -0.11629685 0.08126110 -0.32782773 -0.4397657
## M 0.1389117 -0.08458265 -0.14565003 -0.72512857 -0.2981872
## Sim -0.1500914 0.06421619 0.10816066 0.53220067 0.3842927
## Sio 0.0443487 -0.20414466 0.01018199 -0.27022990 -0.1986540
## V -0.2351446 -0.03694406 -0.19258292 0.07425101 -0.2399229
## ZRo 0.1406930 -0.46339551 0.05857676 -0.05581753 -0.4161484
## ZWi -1.1282488 0.33360379 -0.94134901 0.64348154 0.4876927
#bet2$lw : the weights of the groups are the relative frequencies of the categorical variable
#summary(betw$site)/length(betw$site)
bet2$ratio #the between group variance is the inertia of the between PCA
## [1] 0.1430046
#The between inertia is here to 0.1430 i.e. 14.3% of the total inertia is due to the factor.
plot(bet2)
# VIF coefficients with "corvif" function (Zuur et al. 2009, Highland Statistics LTD)
corvif(vif)
##
##
## Variance inflation factors
##
## GVIF
## altitude 26.904671
## sfroy 6.125720
## ddeg 19.933109
## mbal7 19.525335
## pday 11.213173
## srad7 2.289304
#removing "altitude""
vif2<-vif[,-1]
corvif(vif2)
##
##
## Variance inflation factors
##
## GVIF
## sfroy 4.947536
## ddeg 3.455439
## mbal7 19.507059
## pday 11.195943
## srad7 2.288910
#removing "mbal7"
vif3<-vif2[,-3]
corvif(vif3)
##
##
## Variance inflation factors
##
## GVIF
## sfroy 2.503023
## ddeg 2.641762
## pday 1.297041
## srad7 1.183373
#all are now <3
#PCA on both matrices
dudi.plant<-dudi.pca(plant,scale=TRUE,scan=FALSE,nf=2)
dudi.env<-dudi.pca(env,scale=TRUE,scan=FALSE,nf=2)
#Relative variation of EV
dudi.plant$eig/sum(dudi.plant$eig)
## [1] 0.585755000 0.342009461 0.041855924 0.020398002 0.009981613
dudi.env$eig/sum(dudi.env$eig)
## [1] 0.617638103 0.195731103 0.135461185 0.042516027 0.005120133 0.003533449
#Equal row weights in the 2 analyses ?
all.equal(dudi.plant$lw,dudi.env$lw)
## [1] TRUE
#Must be equal in the two separate ordinations
#Co-inertia analysis
coia.plant.env<-coinertia(dudi.plant,dudi.env,scan=FALSE,nf=2)
#Relative variation on 1st EV
coia.plant.env$eig[1]/sum(coia.plant.env$eig)
## [1] 0.9538968
summary(coia.plant.env)
## Coinertia analysis
##
## Class: coinertia dudi
## Call: coinertia(dudiX = dudi.plant, dudiY = dudi.env, scannf = FALSE,
## nf = 2)
##
## Total inertia: 5.891
##
## Eigenvalues:
## Ax1 Ax2 Ax3 Ax4 Ax5
## 5.620e+00 2.651e-01 6.175e-03 3.127e-04 1.475e-05
##
## Projected inertia (%):
## Ax1 Ax2 Ax3 Ax4 Ax5
## 9.539e+01 4.500e+00 1.048e-01 5.308e-03 2.503e-04
##
## Cumulative projected inertia (%):
## Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
## 95.39 99.89 99.99 100.00 100.00
##
## Eigenvalues decomposition:
## eig covar sdX sdY corr
## 1 5.6198688 2.370626 1.642021 1.920674 0.7516759
## 2 0.2651137 0.514892 1.323993 1.046312 0.3716799
##
## Inertia & coinertia X (dudi.plant):
## inertia max ratio
## 1 2.696234 2.928775 0.9206013
## 12 4.449191 4.638822 0.9591209
##
## Inertia & coinertia Y (dudi.env):
## inertia max ratio
## 1 3.688989 3.705829 0.9954559
## 12 4.783759 4.880215 0.9802352
##
## RV:
## 0.4354129
#MC procedure to test whether the co-inertia between the two multivariate data sets significantly deviate from chance
randtest(coia.plant.env,nrepet=999)
## Monte-Carlo test
## Call: randtest.coinertia(xtest = coia.plant.env, nrepet = 999)
##
## Observation: 0.4354129
##
## Based on 999 replicates
## Simulated p-value: 0.001
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.594672e+01 4.010361e-02 6.145131e-04
plot(coia.plant.env)
tab.mfa <- MFA(tab, group=grn, type=c("c","c","n"), ncp = 2, name.group = c("Env", "Plant", "Site"), num.group.sup = 3)
tab.mfa$group
## $Lg
## Env Plant Site MFA
## Env 1.00177882 0.04453477 0.3777814 0.9011226
## Plant 0.04453477 1.01195212 0.1223304 0.9098842
## Site 0.37778139 0.12233036 8.0000000 0.4307141
## MFA 0.90112256 0.90988416 0.4307141 1.5597035
##
## $RV
## Env Plant Site MFA
## Env 1.00000000 0.04423167 0.13344726 0.7209031
## Plant 0.04423167 1.00000000 0.04299414 0.7242443
## Site 0.13344726 0.04299414 1.00000000 0.1219335
## MFA 0.72090313 0.72424432 0.12193350 1.0000000
##
## $coord
## Dim.1 Dim.2
## Env 0.5632296 0.4542213
## Plant 0.5978928 0.4068062
##
## $contrib
## Dim.1 Dim.2
## Env 48.50734 52.7534
## Plant 51.49266 47.2466
##
## $cos2
## Dim.1 Dim.2
## Env 0.3166643 0.2059506
## Plant 0.3532536 0.1635367
##
## $dist2
## [1] 1.001779 1.011952
##
## $correlation
## Dim.1 Dim.2
## Env 0.7554947 0.6811982
## Plant 0.7744394 0.6401887
##
## $coord.sup
## Dim.1 Dim.2
## Site 0.1916033 0.2849254
##
## $contrib.sup
## Dim.1 Dim.2
## Site 16.50156 33.09132
##
## $cos2.sup
## Dim.1 Dim.2
## Site 0.004588976 0.01014781
##
## $dist2.sup
## [1] 8
tab.mfa$quanti.var
## $coord
## Dim.1 Dim.2
## sfroy -1.049058e+02 -6.586679e+01
## ddeg 1.981954e+02 1.973658e+02
## mbal7 -1.308467e+01 -1.200325e+01
## pday -1.903840e+00 -3.602819e-01
## srad7 -1.696148e+03 1.522292e+03
## nF -9.586408e-01 3.308713e-02
## nB 6.249901e+00 5.308959e+00
## lF -6.706636e-02 2.312110e-04
## lB 8.205992e-01 4.493109e-01
## h 3.523836e+00 2.779313e+00
##
## $contrib
## Dim.1 Dim.2
## sfroy 1.823590e-01 9.694433e-02
## ddeg 6.509017e-01 8.704284e-01
## mbal7 2.836963e-03 3.219492e-03
## pday 6.006057e-05 2.900510e-06
## srad7 4.767119e+01 5.178281e+01
## nF 8.915845e-01 1.432286e-03
## nB 3.789632e+01 3.687487e+01
## lF 4.363754e-03 6.994051e-08
## lB 6.533004e-01 2.641226e-01
## h 1.204709e+01 1.010617e+01
##
## $cos2
## Dim.1 Dim.2
## sfroy 0.26292176 1.036479e-01
## ddeg 0.20382609 2.021233e-01
## mbal7 0.21541724 1.812813e-01
## pday 0.28109874 1.006660e-02
## srad7 0.55368310 4.459949e-01
## nF 0.12858371 1.531765e-04
## nB 0.56747460 4.094670e-01
## lF 0.07649039 9.091058e-07
## lB 0.60570387 1.815900e-01
## h 0.49959132 3.107838e-01
##
## $cor
## Dim.1 Dim.2
## sfroy -0.5127590 -0.3219439998
## ddeg 0.4514710 0.4495812856
## mbal7 -0.4641306 -0.4257714391
## pday -0.5301875 -0.1003324431
## srad7 -0.7440989 0.6678285151
## nF -0.3585857 0.0123764514
## nB 0.7533091 0.6398961060
## lF -0.2765690 0.0009534704
## lB 0.7782698 0.4261338290
## h 0.7068177 0.5574799038
#group : number of variables in each group, type : c for quanti. variables, n for categorical ones, ncp : number of dimensions, num.group.sup : index of the illustrative groupe.
plot(tab.mfa, choix="ind", habillage=11)
plotellipses(tab.mfa, keepvar=11)
coeffRV(plant,env)
## $rv
## [1] 0.04423167
##
## $rvstd
## [1] 0.934212
##
## $mean
## [1] 0.02129687
##
## $variance
## [1] 0.0006026968
##
## $skewness
## [1] 2.57552
##
## $p.value
## [1] 0.1327553
Escofier & Pagès (2008) * Ellipse de confiance pour traduire la variabilité des individus autour des centres de gravité, analogue bi-dimensionnel de l’intervalle de confiance que l’on calcule usuellement autour d’une moyenne. * Procédure “bootstrap” (pour chaque échantillon “bootstrap”, calculer les centres de gravité et les projeter en supplémentaire sur les plans). L’ellipse est centrée sur le centre de gravité initial et contient 95 % des n centres de gravité bootstrap. * La taille d’une ellipse ainsi obtenue dépend de la variabilité (dans le plan factoriel) des individus présentant la modalité étudiée mais aussi de son effectif.
tab.mfa <- MFA(tab, group=grn, type=c("c","c","n"), ncp = 2, name.group = c("Env", "Plant", "Site"), num.group.sup = 3)
tab.mfa$group
## $Lg
## Env Plant Site MFA
## Env 1.00578535 0.06079793 0.3859738 0.9111682
## Plant 0.06079793 1.01195212 0.1223304 0.9164364
## Site 0.38597378 0.12233036 8.0000000 0.4342376
## MFA 0.91116824 0.91643643 0.4342376 1.5612989
##
## $RV
## Env Plant Site MFA
## Env 1.00000000 0.06026375 0.13606930 0.7271146
## Plant 0.06026375 1.00000000 0.04299414 0.7290870
## Site 0.13606930 0.04299414 1.00000000 0.1228682
## MFA 0.72711464 0.72908696 0.12286818 1.0000000
##
## $coord
## Dim.1 Dim.2
## Env 0.5460787 0.4890057
## Plant 0.6244881 0.3791589
##
## $contrib
## Dim.1 Dim.2
## Env 46.65079 56.32638
## Plant 53.34921 43.67362
##
## $cos2
## Dim.1 Dim.2
## Env 0.2964866 0.2377511
## Plant 0.3853793 0.1420635
##
## $dist2
## [1] 1.005785 1.011952
##
## $correlation
## Dim.1 Dim.2
## Env 0.7493795 0.7123192
## Plant 0.7913168 0.6184104
##
## $coord.sup
## Dim.1 Dim.2
## Site 0.1830101 0.2957206
##
## $contrib.sup
## Dim.1 Dim.2
## Site 15.63431 34.06273
##
## $cos2.sup
## Dim.1 Dim.2
## Site 0.004186587 0.01093133
##
## $dist2.sup
## [1] 8
tab.mfa$quanti.var
## $coord
## Dim.1 Dim.2
## altitude -227.7151042 -1.931874e+02
## sfroy -108.1251688 -6.453445e+01
## ddeg 207.5380379 1.962453e+02
## mbal7 -13.5659198 -1.180570e+01
## pday -1.9175757 -3.118802e-01
## srad7 -1652.9160553 1.569105e+03
## nF -0.9463669 7.309893e-02
## nB 6.3860249 5.119145e+00
## lF -0.0658663 3.506274e-03
## lB 0.8326324 4.258233e-01
## h 3.6141269 2.695553e+00
##
## $contrib
## Dim.1 Dim.2
## altitude 8.520931e-01 8.269039e-01
## sfroy 1.921131e-01 9.227429e-02
## ddeg 7.077808e-01 8.532886e-01
## mbal7 3.024140e-03 3.088024e-03
## pday 6.042385e-05 2.155129e-06
## srad7 4.489572e+01 5.455082e+01
## nF 8.618896e-01 6.933439e-03
## nB 3.924585e+01 3.400333e+01
## lF 4.175025e-03 1.595212e-05
## lB 6.671741e-01 2.352803e-01
## h 1.257012e+01 9.428062e+00
##
## $cos2
## Dim.1 Dim.2
## altitude 0.28534915 0.2053763818
## sfroy 0.27930664 0.0994972127
## ddeg 0.22349521 0.1998348934
## mbal7 0.23155471 0.1753631728
## pday 0.28516940 0.0075435127
## srad7 0.52581786 0.4738464197
## nF 0.12531217 0.0007476464
## nB 0.59246318 0.3807105626
## lF 0.07377749 0.0002090687
## lB 0.62359816 0.1631011763
## h 0.52552123 0.2923339587
##
## $cor
## Dim.1 Dim.2
## altitude -0.5341808 -0.45318471
## sfroy -0.5284947 -0.31543179
## ddeg 0.4727528 0.44702896
## mbal7 -0.4812013 -0.41876386
## pday -0.5340125 -0.08685340
## srad7 -0.7251330 0.68836503
## nF -0.3539946 0.02734312
## nB 0.7697163 0.61701747
## lF -0.2716201 0.01445921
## lB 0.7896823 0.40385787
## h 0.7249284 0.54067916
#group : number of variables in each group, type : c for quanti. variables, n for categorical ones, ncp : number of dimensions, num.group.sup : index of the illustrative groupe.
plot(tab.mfa, choix="ind", habillage=12)
plotellipses(tab.mfa, keepvar=12)
# The similarity between the geometrical representations derived from each group of variables is measured by the RV coefficient
coeffRV(plant,env)
## $rv
## [1] 0.06026375
##
## $rvstd
## [1] 1.561259
##
## $mean
## [1] 0.02196013
##
## $variance
## [1] 0.0006019073
##
## $skewness
## [1] 2.560513
##
## $p.value
## [1] 0.07485604
(=regression followed by PCA)
Borcard et al. (2011)
Legendre & Legendre (2012)
#For RDA : do not take derived measures
#Also take altitude as environmental variables (to do enventuallay variance partitioning later)
#Standardize and convert into data frame
mes<-as.data.frame(scale(mes))
env<-as.data.frame(scale(env))
# rda(Y,X,W) where Y is the response matrix,X is the matrix of explanatory variables and W is an optional matrix of covariables
rda.site<-rda(mes~.,env) #same as : rda(mes,env), but need formula for anova
#summary(rda.site)
coef(rda.site)
## RDA1 RDA2 RDA3 RDA4 RDA5
## altitude -0.1461124301 -0.212946868 0.1678457 0.04594887 -0.57968534
## sfroy -0.0229693194 0.200755604 0.2049735 -0.05977453 0.14824250
## ddeg -0.0502858736 -0.058323382 0.1613492 0.24538846 -0.48642561
## mbal7 0.0391286463 -0.185232505 -0.2762264 0.45396415 -0.02885589
## pday -0.0841985522 0.208001894 0.1400771 -0.28328422 0.06435053
## srad7 0.0008958228 -0.005891491 -0.1189128 0.05192108 -0.07594276
plot(rda.site)
# percentage of variance explained by axis 1
rda.site$CCA$eig[1]/rda.site$tot.chi
## RDA1
## 0.3215998
# percentage of variance explained by axis 2
rda.site$CCA$eig[2]/rda.site$tot.chi
## RDA2
## 0.07405687
# R-squared and adjusted-R2
(R2 <- RsquareAdj(rda.site)$r.squared)
## [1] 0.4004635
(R2_adj <- RsquareAdj(rda.site)$adj.r.squared)
## [1] 0.3285192
# check Variance Inflation Factors (colinearity between predictors) and simplify the rda-model
forward.sel(mes, env, adjR2thresh=R2_adj)
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.078000 (superior to 0.050000)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 altitude 1 0.26796289 0.2679629 0.2546531 20.132803 0.001
## 2 pday 5 0.08066717 0.3486301 0.3245052 6.687485 0.004
vif.cca(rda.site)
## altitude sfroy ddeg mbal7 pday srad7
## 26.904671 6.125720 19.933109 19.525335 11.213173 2.289304
#analyse linear dependencies among constraints and conditions
#VIF are the inverse of tolerance, > 5, or worse, 10 indicate collinearity
anova(rda.site,by="axis")
## Model: rda(formula = mes ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## RDA1 1 1.60800 27.3571 199 0.005 **
## RDA2 1 0.37028 6.2997 199 0.005 **
## RDA3 1 0.02215 0.3768 99 0.670
## RDA4 1 0.00150 0.0256 99 1.000
## RDA5 1 0.00038 0.0065 99 1.000
## Residual 51 2.99768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda.site,by="term")
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
##
## Model: rda(formula = mes ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## altitude 1 1.33981 22.3475 99 0.01 **
## sfroy 1 0.13144 2.1923 99 0.10 .
## ddeg 1 0.03743 0.6244 99 0.52
## mbal7 1 0.20640 3.4426 99 0.04 *
## pday 1 0.27908 4.6549 99 0.03 *
## srad7 1 0.00816 0.1361 99 0.92
## Residual 50 2.99768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#For RDA : do not take derived measures
#Also take altitude as environmental variables (to do enventuallay variance partitioning later)
#Standardize and convert into data frame
mes<-as.data.frame(scale(mes))
env<-as.data.frame(scale(env))
# rda(Y,X,W) where Y is the response matrix,X is the matrix of explanatory variables and W is an optional matrix of covariables
rda.site<-rda(mes~.,env) #same as : rda(mes,env), but need formula for anova
#summary(rda.site)
coef(rda.site)
## RDA1 RDA2 RDA3 RDA4 RDA5
## sfroy -0.0576900845 0.160945595 0.23726309 -0.04810805 -0.02256243
## ddeg 0.0632530523 0.120033745 0.01645905 0.20662753 0.01814911
## mbal7 0.0397618664 -0.199143491 -0.26643542 0.45378109 0.17377021
## pday -0.0937937651 0.206098043 0.13637959 -0.28202431 -0.22475938
## srad7 0.0003965869 -0.007692333 -0.12219933 0.04889505 0.15326661
plot(rda.site)
# percentage of variance explained by axis 1
rda.site$CCA$eig[1]/rda.site$tot.chi
## RDA1
## 0.3077382
# percentage of variance explained by axis 2
rda.site$CCA$eig[2]/rda.site$tot.chi
## RDA2
## 0.06669089
# R-squared and adjusted-R2
(R2 <- RsquareAdj(rda.site)$r.squared) #compare with the R^2 obtained with Mantel test
## [1] 0.3788682
(R2_adj <- RsquareAdj(rda.site)$adj.r.squared)
## [1] 0.317973
# check Variance Inflation Factors (colinearity between predictors) and simplify the rda-model
library(packfor)
forward.sel(mes, env, adjR2thresh=R2_adj)
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.052000 (superior to 0.050000)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 mbal7 3 0.25176620 0.2517662 0.2381620 18.506437 0.001
## 2 ddeg 2 0.05077202 0.3025382 0.2767063 3.930953 0.018
vif.cca(rda.site)
## sfroy ddeg mbal7 pday srad7
## 4.947536 3.455439 19.507059 11.195943 2.288910
#analyse linear dependencies among constraints and conditions
#VIF are the inverse of tolerance, > 5, or worse, 10 indicate collinearity
anova(rda.site,by="axis")
## Model: rda(formula = mes ~ sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## RDA1 1 1.53869 25.2678 199 0.005 **
## RDA2 1 0.33345 5.4759 199 0.005 **
## RDA3 1 0.02065 0.3392 99 0.770
## RDA4 1 0.00150 0.0246 99 1.000
## RDA5 1 0.00005 0.0008 99 1.000
## Residual 51 3.10566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda.site,by="term")
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
##
## Model: rda(formula = mes ~ sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## sfroy 1 0.97300 15.9783 99 0.01 **
## ddeg 1 0.37843 6.2145 99 0.01 **
## mbal7 1 0.25120 4.1251 99 0.01 **
## pday 1 0.28356 4.6565 99 0.03 *
## srad7 1 0.00815 0.1338 99 0.92
## Residual 51 3.10566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plant <- env_mes[,4:8] #mesures (without "derived"" measures)
env <- env_mes[,10:14] #environmental variables (NB : without altitude)
spa<-spa[,-1] #coordonnées des "pop" (x/y)
#Compute distance matrices to run Mantel test
plant.eu <- vegdist(scale(plant), "euclidean") #computes dissimilarity matrices/indices
env.eu <- vegdist(scale(env), "euclidean")
spa.eu<-vegdist(scale(spa),"euclidean")
#H0 : the 2 matrices are uncorrelated
#mantel(xids,ydis)
#mantel site betw. plant measures & env.variables
(mantel_site<-mantel(plant.eu, env.eu, strata=env_mes$site, method="kendall",permutations=1000))
##
## Mantel statistic based on Kendall's rank correlation tau
##
## Call:
## mantel(xdis = plant.eu, ydis = env.eu, method = "kendall", permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: 0.243
## Significance: 0.000999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0885 0.1029 0.1170 0.1346
##
## Based on 1000 permutations, stratified within env_mes$site
#mantel site betw. plant measures & geographic coordinates
(mantel_site<-mantel(plant.eu, spa.eu, strata=env_mes$site, method="kendall",permutations=1000))
##
## Mantel statistic based on Kendall's rank correlation tau
##
## Call:
## mantel(xdis = plant.eu, ydis = spa.eu, method = "kendall", permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: -0.005196
## Significance: 0.67732
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.00309 0.00443 0.00602 0.00678
##
## Based on 1000 permutations, stratified within env_mes$site
#partial Mantel test :
#Partial Mantel statistic uses partial correlation conditioned on the third matrix. Only the first matrix is permuted so that the correlation structure between second and first matrices is kept constant. Although mantel.partial silently accepts other methods than "pearson", partial correlations will probably be wrong with other methods.
#mantel.partial(xdis, ydis, zdis)
(partial_mantel_site<-mantel.partial(plant.eu,env.eu,spa.eu,strata=env_mes$site,permutations=1000))
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = plant.eu, ydis = env.eu, zdis = spa.eu, permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: 0.4535
## Significance: 0.000999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.229 0.260 0.282 0.297
##
## Based on 1000 permutations, stratified within env_mes$site
Lecture Notes (Kienast et al. 2012)
plot(tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)); text(tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all))
model_alt<-tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(prune.tree(model_alt));
#Additional deviance (fit) achieved by adding more nodes beyond (3)4 seems marginal (cost-complexity curve begins to decline at this point). This suggests that the tree could potentially be pruned to just (3)4 terminal branches without a great loss of predictive power to achieve a more genuine predictive model (Logan 2010).
plot(prune.tree(model_alt,best=4));text(prune.tree(model_alt,best=4))
for.fit_alt <- randomForest(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(for.fit_alt)
varImpPlot(for.fit_alt)
print(for.fit_alt) # view results
##
## Call:
## randomForest(formula = ratiopond ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = all)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 2
##
## Mean of squared residuals: 0.005682513
## % Var explained: 23.84
importance(for.fit_alt) # importance of each predictor
## IncNodePurity
## altitude 1.0168356
## sfroy 0.9476718
## ddeg 0.6329066
## mbal7 0.5333658
## pday 0.3360027
## srad7 0.4900198
plot(tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)); text(tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all))
model<-tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(prune.tree(model));
#The tree could potentially be pruned to just (3)4 terminal branches without a great loss of predictive power to achieve a more genuine predictive model (Logan 2010).
plot(prune.tree(model,best=4));text(prune.tree(model,best=4))
for.fit <- randomForest(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(for.fit)
varImpPlot(for.fit)
print(for.fit) # view results
##
## Call:
## randomForest(formula = ratiopond ~ sfroy + ddeg + mbal7 + pday + srad7, data = all)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 1
##
## Mean of squared residuals: 0.005681186
## % Var explained: 23.86
importance(for.fit) # importance of each predictor
## IncNodePurity
## sfroy 1.0598047
## ddeg 0.8008486
## mbal7 0.7598765
## pday 0.5966725
## srad7 0.7326774
#summary(A)
#hist(A$ratiopond.tr)
plot(A$ratiopond~A$altitude)
bwplot(A$ratiopond~as.factor(A$altitude))
#levels(as.factor(A$altitude))
lm.A <- aov((A$ratiopond.tr)~as.factor(A$altitude))
#summary.lm(lm.A)
TukeyHSD(lm.A)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (A$ratiopond.tr) ~ as.factor(A$altitude))
##
## $`as.factor(A$altitude)`
## diff lwr upr p adj
## 1778-1547 0.012396023 -0.05177806 0.076570108 0.9935819
## 1894-1547 -0.038236302 -0.10241039 0.025937783 0.5224629
## 2001-1547 -0.009993451 -0.07416754 0.054180633 0.9976775
## 2303-1547 -0.005466761 -0.06964085 0.058707324 0.9998768
## 2460-1547 -0.066990020 -0.13116410 -0.002815936 0.0351762
## 1894-1778 -0.050632325 -0.11480641 0.013541759 0.2105483
## 2001-1778 -0.022389474 -0.08656356 0.041784610 0.9156417
## 2303-1778 -0.017862784 -0.08203687 0.046311300 0.9668091
## 2460-1778 -0.079386044 -0.14356013 -0.015211959 0.0061684
## 2001-1894 0.028242851 -0.03593123 0.092416935 0.8018544
## 2303-1894 0.032769541 -0.03140454 0.096943625 0.6828376
## 2460-1894 -0.028753719 -0.09292780 0.035420366 0.7895479
## 2303-2001 0.004526690 -0.05964739 0.068700774 0.9999515
## 2460-2001 -0.056996570 -0.12117065 0.007177515 0.1131046
## 2460-2303 -0.061523260 -0.12569734 0.002650825 0.0685233
qqPlot(resid(lm.A))
plot(aov((A$ratiopond.tr)~as.factor(A$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=A), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=A))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.21679 0.00519 41.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.334 3.766 3.426 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0528 Deviance explained = 7.04%
## GCV = 0.0049673 Scale est. = 0.0048477 n = 180
#k is the dimension of the basis used to represent the smooth term
#if edf=1 => straight line (linear)
#summary(Lie)
#hist(Lie$ratiopond.tr)
plot(Lie$ratiopond~Lie$altitude)
bwplot(Lie$ratiopond~as.factor(Lie$altitude))
#levels(as.factor(Lie$altitude))
lm.Lie <- aov((Lie$ratiopond.tr)~as.factor(Lie$altitude))
#summary.lm(lm.Lie)
TukeyHSD(lm.Lie)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Lie$ratiopond.tr) ~ as.factor(Lie$altitude))
##
## $`as.factor(Lie$altitude)`
## diff lwr upr p adj
## 1602-1088 -0.005191541 -0.070385091 0.06000201 0.9999849
## 1822-1088 0.083225999 0.018032449 0.14841955 0.0035384
## 2151-1088 0.086624529 0.021430979 0.15181808 0.0019950
## 2299-1088 0.049322319 -0.015871231 0.11451587 0.2722959
## 2425-1088 0.073345584 0.008152034 0.13853913 0.0164586
## 2539-1088 0.162001416 0.096807866 0.22719497 0.0000000
## 1822-1602 0.088417540 0.023223990 0.15361109 0.0014618
## 2151-1602 0.091816070 0.026622520 0.15700962 0.0007980
## 2299-1602 0.054513860 -0.010679690 0.11970741 0.1683335
## 2425-1602 0.078537125 0.013343575 0.14373067 0.0075219
## 2539-1602 0.167192957 0.101999407 0.23238651 0.0000000
## 2151-1822 0.003398530 -0.061795020 0.06859208 0.9999988
## 2299-1822 -0.033903680 -0.099097230 0.03128987 0.7148702
## 2425-1822 -0.009880415 -0.075073965 0.05531314 0.9993501
## 2539-1822 0.078775417 0.013581867 0.14396897 0.0072467
## 2299-2151 -0.037302210 -0.102495760 0.02789134 0.6141665
## 2425-2151 -0.013278945 -0.078472495 0.05191461 0.9965568
## 2539-2151 0.075376887 0.010183337 0.14057044 0.0121968
## 2425-2299 0.024023265 -0.041170285 0.08921681 0.9283041
## 2539-2299 0.112679097 0.047485547 0.17787265 0.0000128
## 2539-2425 0.088655832 0.023462282 0.15384938 0.0014020
qqPlot(resid(lm.Lie))
plot(aov((Lie$ratiopond.tr)~as.factor(Lie$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Lie), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Lie))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.229500 0.004855 47.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.942 3.998 20.21 1.69e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.273 Deviance explained = 28.7%
## GCV = 0.0050683 Scale est. = 0.0049491 n = 210
#summary(Lou)
#hist(Lou$ratiopond.tr)
plot(Lou$ratiopond~Lou$altitude)
bwplot(Lou$ratiopond~as.factor(Lou$altitude))
#levels(as.factor(Lou$altitude))
lm.Lou <- aov((Lou$ratiopond.tr)~as.factor(Lou$altitude))
#summary.lm(lm.Lou)
TukeyHSD(lm.Lou)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Lou$ratiopond.tr) ~ as.factor(Lou$altitude))
##
## $`as.factor(Lou$altitude)`
## diff lwr upr p adj
## 2132-1978 -0.0038852560 -0.06563616 0.05786565 0.9997933
## 2308-1978 0.0001450315 -0.06160587 0.06189594 1.0000000
## 2440-1978 -0.0512232579 -0.11297416 0.01052765 0.1535574
## 2590-1978 0.0434273563 -0.01885361 0.10570832 0.3081723
## 2308-2132 0.0040302876 -0.05772062 0.06578119 0.9997610
## 2440-2132 -0.0473380019 -0.10908891 0.01441290 0.2181861
## 2590-2132 0.0473126124 -0.01496835 0.10959358 0.2263027
## 2440-2308 -0.0513682895 -0.11311919 0.01038262 0.1514541
## 2590-2308 0.0432823248 -0.01899864 0.10556329 0.3115261
## 2590-2440 0.0946506143 0.03236965 0.15693158 0.0004451
qqPlot(resid(lm.Lou))
plot(aov((Lou$ratiopond.tr)~as.factor(Lou$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Lou), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Lou))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.248987 0.006045 41.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.728 3.955 4.224 0.00304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0813 Deviance explained = 10.4%
## GCV = 0.0056239 Scale est. = 0.0054455 n = 149
#summary(M)
#hist(M$ratiopond.tr)
plot(M$ratiopond~M$altitude)
bwplot(M$ratiopond~as.factor(M$altitude))
#levels(as.factor(M$altitude))
lm.M <- aov((M$ratiopond.tr)~as.factor(M$altitude))
#summary.lm(lm.M)
TukeyHSD(lm.M)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (M$ratiopond.tr) ~ as.factor(M$altitude))
##
## $`as.factor(M$altitude)`
## diff lwr upr p adj
## 1757-1660 0.014883268 -0.050672339 0.080438875 0.9937640
## 1995-1660 -0.077996858 -0.143552465 -0.012441251 0.0087453
## 2130-1660 -0.055917340 -0.121472947 0.009638267 0.1506685
## 2265-1660 -0.032858475 -0.098414082 0.032697132 0.7489194
## 2420-1660 0.063324322 -0.002231285 0.128879929 0.0658770
## 2520-1660 -0.049565578 -0.115121185 0.015990028 0.2730049
## 1995-1757 -0.092880126 -0.158435733 -0.027324519 0.0007220
## 2130-1757 -0.070800608 -0.136356215 -0.005245001 0.0249847
## 2265-1757 -0.047741743 -0.113297350 0.017813864 0.3171567
## 2420-1757 0.048441054 -0.017114553 0.113996661 0.2997772
## 2520-1757 -0.064448846 -0.130004453 0.001106760 0.0574142
## 2130-1995 0.022079518 -0.043476089 0.087635125 0.9528477
## 2265-1995 0.045138383 -0.020417224 0.110693990 0.3863651
## 2420-1995 0.141321180 0.075765573 0.206876787 0.0000000
## 2520-1995 0.028431279 -0.037124327 0.093986886 0.8551828
## 2265-2130 0.023058865 -0.042496742 0.088614472 0.9421343
## 2420-2130 0.119241662 0.053686055 0.184797269 0.0000035
## 2520-2130 0.006351761 -0.059203845 0.071907368 0.9999520
## 2420-2265 0.096182797 0.030627190 0.161738404 0.0003934
## 2520-2265 -0.016707103 -0.082262710 0.048848504 0.9884316
## 2520-2420 -0.112889900 -0.178445507 -0.047334294 0.0000140
qqPlot(resid(lm.M))
plot(aov((M$ratiopond.tr)~as.factor(M$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=M), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=M))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.208629 0.005056 41.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.912 3.995 9.546 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.15 Deviance explained = 16.6%
## GCV = 0.0054979 Scale est. = 0.0053693 n = 210
#summary(Sim)
#hist(Sim$ratiopond.tr)
plot(Sim$ratiopond~Sim$altitude)
bwplot(Sim$ratiopond~as.factor(Sim$altitude))
#levels(as.factor(Sim$altitude))
lm.Sim <- aov((Sim$ratiopond.tr)~as.factor(Sim$altitude))
#summary.lm(lm.Sim)
TukeyHSD(lm.Sim)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Sim$ratiopond.tr) ~ as.factor(Sim$altitude))
##
## $`as.factor(Sim$altitude)`
## diff lwr upr p adj
## 2089-1977 -0.004816182 -0.05793337 0.048301008 0.9953286
## 2193-1977 0.016660301 -0.03645689 0.069777491 0.8460541
## 2292-1977 -0.051030024 -0.10414721 0.002087166 0.0645029
## 2193-2089 0.021476483 -0.03164071 0.074593673 0.7180608
## 2292-2089 -0.046213842 -0.09933103 0.006903348 0.1116275
## 2292-2193 -0.067690325 -0.12080751 -0.014573135 0.0064849
qqPlot(resid(lm.Sim))
plot(aov((Sim$ratiopond.tr)~as.factor(Sim$altitude)))
plot(gam(ratiopond~s(altitude, k=4), data=Sim), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=4), data=Sim))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.239959 0.006171 38.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 2.526 2.834 3.916 0.0123 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0664 Deviance explained = 8.62%
## GCV = 0.0047081 Scale est. = 0.0045697 n = 120
#summary(Sio)
#hist(Sio$ratiopond.tr)
plot(Sio$ratiopond~Sio$altitude)
bwplot(Sio$ratiopond~as.factor(Sio$altitude))
#levels(as.factor(Sio$altitude))
lm.Sio <- aov((Sio$ratiopond.tr)~as.factor(Sio$altitude))
#summary.lm(lm.Sio)
TukeyHSD(lm.Sio)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Sio$ratiopond.tr) ~ as.factor(Sio$altitude))
##
## $`as.factor(Sio$altitude)`
## diff lwr upr p adj
## 1541-1337 -0.055579036 -0.134452820 0.02329475 0.4058692
## 1600-1337 0.004241479 -0.070961723 0.07944468 1.0000000
## 1777-1337 -0.011883273 -0.087086476 0.06331993 0.9999061
## 2229-1337 0.066886229 -0.008316974 0.14208943 0.1256973
## 2353-1337 0.007458761 -0.067744442 0.08266196 0.9999974
## 2572-1337 0.089128004 0.013924802 0.16433121 0.0077608
## 2674-1337 0.136958609 0.061755407 0.21216181 0.0000012
## 2761-1337 0.120793331 0.045590129 0.19599653 0.0000330
## 1600-1541 0.059820515 -0.019053269 0.13869430 0.3035438
## 1777-1541 0.043695763 -0.035178021 0.12256955 0.7254661
## 2229-1541 0.122465265 0.043591481 0.20133905 0.0000718
## 2353-1541 0.063037797 -0.015835988 0.14191158 0.2367039
## 2572-1541 0.144707040 0.065833256 0.22358082 0.0000010
## 2674-1541 0.192537645 0.113663861 0.27141143 0.0000000
## 2761-1541 0.176372367 0.097498583 0.25524615 0.0000000
## 1777-1600 -0.016124752 -0.091327955 0.05907845 0.9990880
## 2229-1600 0.062644750 -0.012558453 0.13784795 0.1892068
## 2353-1600 0.003217282 -0.071985921 0.07842048 1.0000000
## 2572-1600 0.084886525 0.009683323 0.16008973 0.0142379
## 2674-1600 0.132717130 0.057513928 0.20792033 0.0000029
## 2761-1600 0.116551852 0.041348650 0.19175505 0.0000748
## 2229-1777 0.078769502 0.003566300 0.15397270 0.0321749
## 2353-1777 0.019342034 -0.055861169 0.09454524 0.9966667
## 2572-1777 0.101011277 0.025808075 0.17621448 0.0011999
## 2674-1777 0.148841883 0.073638680 0.22404509 0.0000001
## 2761-1777 0.132676604 0.057473402 0.20787981 0.0000030
## 2353-2229 -0.059427468 -0.134630671 0.01577573 0.2505708
## 2572-2229 0.022241775 -0.052961427 0.09744498 0.9913663
## 2674-2229 0.070072381 -0.005130822 0.14527558 0.0899209
## 2761-2229 0.053907102 -0.021296100 0.12911030 0.3815014
## 2572-2353 0.081669243 0.006466041 0.15687245 0.0220604
## 2674-2353 0.129499849 0.054296646 0.20470305 0.0000057
## 2761-2353 0.113334570 0.038131368 0.18853777 0.0001367
## 2674-2572 0.047830605 -0.027372597 0.12303381 0.5519136
## 2761-2572 0.031665327 -0.043537875 0.10686853 0.9254739
## 2761-2674 -0.016165278 -0.091368481 0.05903792 0.9990712
qqPlot(resid(lm.Sio))
plot(aov((Sio$ratiopond.tr)~as.factor(Sio$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Sio), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Sio))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.248157 0.005036 49.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 2.146 2.619 34.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.252 Deviance explained = 25.8%
## GCV = 0.0068019 Scale est. = 0.0067211 n = 265
#summary(V)
#hist(V$ratiopond.tr)
plot(V$ratiopond~V$altitude)
bwplot(V$ratiopond~as.factor(V$altitude))
#levels(as.factor(V$altitude))
lm.V <- aov((V$ratiopond.tr)~as.factor(V$altitude))
#summary.lm(lm.V)
TukeyHSD(lm.V)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (V$ratiopond.tr) ~ as.factor(V$altitude))
##
## $`as.factor(V$altitude)`
## diff lwr upr p adj
## 1970-1701 0.072008294 0.002268433 0.14174815 0.0379985
## 2192-1701 -0.103997637 -0.173737497 -0.03425778 0.0002926
## 2327-1701 -0.030648891 -0.100388752 0.03909097 0.8472306
## 2429-1701 0.086619673 0.016879812 0.15635953 0.0051112
## 2532-1701 0.099915096 0.030175235 0.16965496 0.0005975
## 2645-1701 0.078825222 0.009085361 0.14856508 0.0156606
## 2192-1970 -0.176005930 -0.245745791 -0.10626607 0.0000000
## 2327-1970 -0.102657185 -0.172397046 -0.03291732 0.0003709
## 2429-1970 0.014611379 -0.055128482 0.08435124 0.9959774
## 2532-1970 0.027906802 -0.041833059 0.09764666 0.8966579
## 2645-1970 0.006816928 -0.062922933 0.07655679 0.9999494
## 2327-2192 0.073348745 0.003608884 0.14308861 0.0321510
## 2429-2192 0.190617310 0.120877449 0.26035717 0.0000000
## 2532-2192 0.203912732 0.134172871 0.27365259 0.0000000
## 2645-2192 0.182822858 0.113082997 0.25256272 0.0000000
## 2429-2327 0.117268564 0.047528703 0.18700843 0.0000244
## 2532-2327 0.130563987 0.060824126 0.20030385 0.0000016
## 2645-2327 0.109474113 0.039734252 0.17921397 0.0001082
## 2532-2429 0.013295423 -0.056444438 0.08303528 0.9976160
## 2645-2429 -0.007794451 -0.077534312 0.06194541 0.9998889
## 2645-2532 -0.021089874 -0.090829735 0.04864999 0.9721364
qqPlot(resid(lm.V))
plot(aov((V$ratiopond.tr)~as.factor(V$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=V), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=V))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.267446 0.005389 49.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.977 4 28.86 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.349 Deviance explained = 36.2%
## GCV = 0.0062476 Scale est. = 0.0060995 n = 210
#summary(ZRo)
#hist(ZRo$ratiopond.tr)
plot(ZRo$ratiopond~ZRo$altitude)
bwplot(ZRo$ratiopond~as.factor(ZRo$altitude))
#levels(as.factor(ZRo$altitude))
lm.ZRo <- aov((ZRo$ratiopond.tr)~as.factor(ZRo$altitude))
#summary.lm(lm.ZRo)
TukeyHSD(lm.ZRo)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (ZRo$ratiopond.tr) ~ as.factor(ZRo$altitude))
##
## $`as.factor(ZRo$altitude)`
## diff lwr upr p adj
## 2374-1846 -0.0255357363 -0.10149423 0.050422759 0.9532641
## 2541-1846 0.0230309928 -0.05292750 0.098989488 0.9717680
## 2732-1846 -0.0381892829 -0.11414778 0.037769213 0.7462030
## 2891-1846 -0.0253784975 -0.10133699 0.050579998 0.9546236
## 3071-1846 0.0607720781 -0.01518642 0.136730574 0.2113174
## 3212-1846 -0.1172068476 -0.19316534 -0.041248352 0.0001524
## 2541-2374 0.0485667291 -0.02739177 0.124525225 0.4799366
## 2732-2374 -0.0126535467 -0.08861204 0.063304949 0.9988851
## 2891-2374 0.0001572387 -0.07580126 0.076115734 1.0000000
## 3071-2374 0.0863078143 0.01034932 0.162266310 0.0147919
## 3212-2374 -0.0916711113 -0.16762961 -0.015712616 0.0073565
## 2732-2541 -0.0612202757 -0.13717877 0.014738220 0.2038223
## 2891-2541 -0.0484094903 -0.12436799 0.027549005 0.4840126
## 3071-2541 0.0377410853 -0.03821741 0.113699581 0.7565850
## 3212-2541 -0.1402378403 -0.21619634 -0.064279345 0.0000024
## 2891-2732 0.0128107854 -0.06314771 0.088769281 0.9988041
## 3071-2732 0.0989613610 0.02300287 0.174919857 0.0026591
## 3212-2732 -0.0790175646 -0.15497606 -0.003059069 0.0355393
## 3071-2891 0.0861505756 0.01019208 0.162109071 0.0150878
## 3212-2891 -0.0918283500 -0.16778685 -0.015869854 0.0072026
## 3212-3071 -0.1779789256 -0.25393742 -0.102020430 0.0000000
qqPlot(resid(lm.ZRo))
plot(aov((ZRo$ratiopond.tr)~as.factor(ZRo$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=ZRo), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=ZRo))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.26385 0.00634 41.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.842 3.982 3.876 0.00473 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0628 Deviance explained = 8%
## GCV = 0.0086396 Scale est. = 0.0084404 n = 210
#summary(ZWi)
#hist(ZWi$ratiopond.tr)
plot(ZWi$ratiopond~ZWi$altitude)
bwplot(ZWi$ratiopond~as.factor(ZWi$altitude))
#levels(as.factor(ZWi$altitude))
lm.ZWi <- aov((ZWi$ratiopond.tr)~as.factor(ZWi$altitude))
#summary.lm(lm.ZWi)
TukeyHSD(lm.ZWi)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (ZWi$ratiopond.tr) ~ as.factor(ZWi$altitude))
##
## $`as.factor(ZWi$altitude)`
## diff lwr upr p adj
## 2057-1939 0.022810101 -0.026007424 0.07162763 0.6973310
## 2321-1939 0.006208213 -0.042609312 0.05502574 0.9967018
## 2584-1939 0.057942516 0.009124990 0.10676004 0.0112562
## 2936-1939 0.105389631 0.056572105 0.15420716 0.0000002
## 2321-2057 -0.016601888 -0.065419414 0.03221564 0.8810249
## 2584-2057 0.035132415 -0.013685111 0.08394994 0.2770457
## 2936-2057 0.082579530 0.033762004 0.13139706 0.0000653
## 2584-2321 0.051734303 0.002916777 0.10055183 0.0319289
## 2936-2321 0.099181418 0.050363892 0.14799894 0.0000010
## 2936-2584 0.047447115 -0.001370411 0.09626464 0.0611403
qqPlot(resid(lm.ZWi))
plot(aov((ZWi$ratiopond.tr)~as.factor(ZWi$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=ZWi), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=ZWi))
##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.254376 0.005002 50.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 1.948 2.354 19.09 8.66e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.227 Deviance explained = 23.7%
## GCV = 0.003829 Scale est. = 0.0037537 n = 150
Packages lme4 (Bates et al. 2014) and nlme (Pinheiro et al. 2015)
Theoretical aspects : West et al. (2007)
General equation :
Ri=Xi x \(\beta\) + Zi x bi + \(\epsilon\)i
Because the random effects come from a large population, there is no much point in concentrating on estimating means of our small subset of factor levels […]. Much better to recognize them for what they are, random samples from a much larger population, and to concentrate on their variance. This is the added variation caused by differences between levels of the random effects (Crawley 2007).
Correlation between the pseudoreplicates within a group causes shrinkage. The BLUPs (parameter estimates, ai) are smaller than the fixed effect size (yi--\(\mu\)). When \(\sigma\)a2 (between-group variance which introduces the correlation between the pseudoreplicates within each group) is estimated to be large compared with the estimate of \(\sigma\)2/n (where \(\sigma\) is the residual variance), then the fixed effects and the BLUP are similar. When \(\sigma\)a2 is estimated to be small compared with the estimate of \(\sigma\)2/n, then the fixed effects and the BLUP can be very different. BLUPs given by : ai=(effect size)*(\(\sigma\)a2/(\(\sigma\)a2+\(\sigma\)2/n)) (Crawley 2007).
In the LMM, we estimate the fixed-effect parameters, and the covariance parameters. 2 common methods to estimate these parameters :
Maximum likelihood (ML) estimation : method of obtaining estimates of unknown parameters by optimizing a likelihood function. To apply ML estimation, we first construct the likelihood as a function of the parameters in the specified model, based on distributional assumptions. The maximum likelihood estimates (MLEs) of the parameters are the values of the arguments that maximize the likelihood function (i.e., the values of the parameters that make the observed values of the dependent variable most likely, given the distributional assumptions).
Restricted maximum likelihood (REML) estimation : (sometimes called residual maximum likelihood estimation) was introduced as a method of estimating variance components in the context of unbalanced incomplete block designs. REML is often preferred to ML estimation, because it produces unbiased estimates of covariance parameters by taking into account the loss of degrees of freedom that results from estimating the fixed effects.
To compare models with nested fixed effects but with the same random structure : ML estimation.
In general terms, we use maximum likelihood methods (either REML or ML estimation) to obtain estimates of the covariance parameters in an LMM. We then obtain estimates of the fixed-effect parameters using results from generalized least squares. However, ML estimates of the covariance parameters are biased, whereas REML estimates are not. Note that the variances of the estimated fixed effects are biased downward in both ML and REML estimation.
When used to estimate the covariance parameters, ML and REML estimation are computationally intensive ; both involve the optimization of some objective function, which generally requires starting values for the parameter estimates and several subsequent iterations to find the values of the parameters that maximize the likelihood function.
The information criteria (sometimes referred to as fit criteria) provide a way to assess the fit of a model based on its optimum log-likelihood value, after applying a penalty for the parameters that are estimated in fitting the model. A key feature of the information criteria is that they provide a way to compare any two models fitted to the same set of observations; i.e., the models do not need to be nested. We use the “smaller is better” form for the information criteria discussed in this section; that is, a smaller value of the criterion indicates a “better” fit. Recent work suggests that no one information criterion stands apart as the best criterion to be used when selecting LMMs.
Schwarz criterion (SC) (also known as the Bayesian information criterion - BIC): a model selection criterion designed to find the most probable model (from a Bayesian perspective) given the data. Superficially similar to AICc, SC has 2 components: (i) negative log- likelihood, which measures lack of fit, (ii) a penalty term that varies as a function of sample size and the number of model parameters. SC is equivalent (under certain conditions) to the natural logarithm of the Bayes factor.
In our case, we identify :
#Fit the initial "unconditional" (variance components) model and decide whether to omit the random effects of the level 2 (population).
#site<-all$site
#pop<-all$site
model1<-lmer(ratiopond.tr~altitude+(1|site/pop),all,REML=TRUE)
summary(model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ratiopond.tr ~ altitude + (1 | site/pop)
## Data: all
##
## REML criterion at convergence: -3315.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9186 -0.6229 -0.0130 0.5921 4.2964
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:site (Intercept) 0.002189 0.04678
## site (Intercept) 0.000000 0.00000
## Residual 0.007638 0.08739
## Number of obs: 1704, groups: pop:site, 57; site, 9
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.779e-01 3.436e-02 10.996
## altitude 5.896e-05 1.528e-05 3.857
##
## Correlation of Fixed Effects:
## (Intr)
## altitude -0.982
anova(model1)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## altitude 1 0.11365 0.11365 14.88
#Model validation
plot(model1,type=c("p","smooth")) #fitted vs. residual
plot(model1,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth")) #scale-location
#qqmath(model1,id=0.05) # quantile-quantile
#R2
r.squaredGLMM(model1) #conditional and marginal coefficient of determination for GLMM
## R2m R2c
## 0.0608076 0.2700072
# marginal : concerned with variance explained by fixed factors
# conditional : concerned with variance explained by both random and fixed factors
r.squaredLR(model1) #calculate a coefficient of determination based on the likelihood-ratio test (R_LR²)
## [1] 0.2062804
## attr(,"adj.r.squared")
## [1] -0.04437533
plot(model1,type=c("p","smooth")) #fitted vs. residual
plot(model1,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth")) #scale-location
#qqmath(m1,id=0.05) # quantile-quantile
| Site | Altitude | Collected spikes |
|---|---|---|
| Anzère (A) | ||
| +A1 | 1547 | 30 |
| +A2 | 1778 | 30 |
| +A3 | 1894 | 30 |
| +A4 | 2001 | 30 |
| +A5 | 2303 | 30 |
| +A6 | 2460 | 30 |
| Lienne (Lie) | ||
| +Lie1 | 1088 | 30 |
| +Lie2 | 1602 | 30 |
| +Lie3 | 1822 | 30 |
| +Lie4 | 2151 | 30 |
| +Lie5 | 2299 | 30 |
| +Lie6 | 2425 | 30 |
| +Lie7 | 2539 | 30 |
| Lourtier (Lou) | ||
| +Lou1 | 1978 | 30 |
| +Lou2 | 2132 | 30 |
| +Lou3 | 2308 | 30 |
| +Lou4 | 2440 | 30 |
| +Lou5 | 2590 | 29 |
| Montana (M) | ||
| +M1 | 1660 | 30 |
| +M2 | 1757 | 30 |
| +M3 | 1995 | 30 |
| +M4 | 2130 | 30 |
| +M5 | 2265 | 30 |
| +M6 | 2420 | 30 |
| +M7 | 2520 | 30 |
| Simplon (Sim) | ||
| +Sim1 | 1977 | 30 |
| +Sim2 | 2089 | 30 |
| +Sim3 | 2193 | 30 |
| +Sim4 | 2292 | 30 |
| Sionne (Sio) | ||
| +Sio1 | 1337 | 30 |
| +Sio2 | 1541 | 25 |
| +Sio3 | 1600 | 30 |
| +Sio4 | 1777 | 30 |
| +Sio5 | 2229 | 30 |
| +Sio6 | 2353 | 30 |
| +Sio7 | 2572 | 30 |
| +Sio8 | 2674 | 30 |
| +Sio9 | 2761 | 30 |
| Vercorin (V) | ||
| +V1 | 1701 | 30 |
| +V2 | 1970 | 30 |
| +V3 | 2192 | 30 |
| +V4 | 2327 | 30 |
| +V5 | 2429 | 30 |
| +V6 | 2532 | 30 |
| +V7 | 2645 | 30 |
| Zermatt-Rothorn (ZRo) | ||
| +ZRo1 | 1846 | 30 |
| +ZRo2 | 2374 | 30 |
| +ZRo3 | 2541 | 30 |
| +ZRo4 | 2732 | 30 |
| +ZRo5 | 2891 | 30 |
| +ZRo6 | 3071 | 30 |
| +ZRo7 | 3212 | 30 |
| Zermatt-Wisshorn (ZWi) | ||
| +ZWi1 | 1939 | 30 |
| +ZWi2 | 2057 | 30 |
| +ZWi3 | 2321 | 30 |
| +ZWi4 | 2584 | 30 |
| +ZWi5 | 2936 | 30 |
| ———————- | ——— | —————– |
## Loading required package: RgoogleMaps